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Abstract 

The analysis of space-time data from complex, real-life phenomena re¬ 
quires the use of flexible and physically motivated covariance functions. In 
most cases, it is not possible to explicitly solve the equations of motion for 
the fields or the respective covariance functions. In the statistical litera¬ 
ture, covariance functions are often based on mathematical constructions. 

We propose deriving space-time covariance functions by solving “effective 
equations of motion”, which can be used as statistical representations of 
systems with diffusive behavior. In particular, we propose using the linear 
response theory to formulate space-time covariance functions based on an 
equilibrium effective Hamiltonian. The effective space-time dynamics are 
then generated by a stochastic perturbation around the equilibrium point 
of the classical field Hamiltonian leading to an associated Langevin equa¬ 
tion. We employ a Hamiltonian which extends the classical Gaussian held 
theory by including a curvature term and leads to a diffusive Langevin 
equation. Finally, we derive new forms of space-time covariance functions. 

Keywords: fluctuation dissipation theorem; interpolation; simulation; space- 
time system; maximum entropy. 


1 Introduction 

Covariance functions play a significant role in the analysis of space-time data 
with geostatistical and machine learning methods mm w, in inverse mod¬ 
elling [75], and in data assimilation [21]. Thus, there is active interest in ap¬ 
plications of space-time data analysis and the development of new covariance 
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models mm • Space-time covariance functions commonly used are straight¬ 
forward extensions of purely spatial models (e.g., exponential, Gaussian) and 
constructions based on linear mixtures [39j . More complex models are generated 
using mathematical arguments, e.g. permissibility conditions 140] and simpli¬ 
fying assumptions such as separability of space and time components without 
clear physical motivation for the functional form or the parameters. Several au¬ 
thors have proposed more realistic, non-separable covariance models B3ES3 El- 
Simple space-time covariance models are determined by the variance and con¬ 
stant correlation length and correlation time. In purely spatial cases, richer 
parametric families such as the Whittle-Matern covariance model [H3 ED] or 
the Spartan covariance functions [13 El EH offer more flexibility. The former 
is derived from a fractional stochastic partial differential equation driven by 
white noise [53] EE] • The latter is obtained from a Gaussian field theory with a 
curvature term which also leads to a fourth-order stochastic partial differential 
equation mm- 

There is a long history of ideas from physics that find applications in in¬ 
formation theory and data analysis, including the maximum entropy formalism 
developed by Jaynes [32], more recently the Bayesian field theory 133, and 
machine learning applications of the variational approximation [22]. The con¬ 
nection between statistical physics and space-time statistics on the other hand 
has not been duly appreciated, albeit the former can provide useful models for 
the analysis of space-time data. In statistical physics the literature on space- 
time fields and correlation functions is extensive P3I- For example, the Gaussian 
model of classical field theory is equivalent to a Gaussian random field with a 
specific covariance function pH [43]. 

A formal difference between random fields in statistical field theories and 
those in spatial statistics is that the former are defined by means of local inter¬ 
actions, whereas the latter are defined by means of covariance matrices 0 ID]. 
Both approaches are formally equivalent provided that the covariance function 
of the local interaction model can be explicitly expressed. The interaction- 
based formalism has several advantages for parameter estimation, interpolation 
and simulation which derive from the sparseness of the inverse covariance ma¬ 
trix in the local interaction representation EH ED]- Spartan Spatial Random 
Fields (SSRFs) are based on the local interaction framework in the static (time- 
independent) case. SSRFs admit explicit relations for the covariance function 
in one, two, and three dimensions EH El QUEUES- A similar framework was 
independently proposed by Farmer m- More recently, covariance functions 
similar to SSRF were derived for meteorological applications from polynomials 
of the diffusion operator EE:, ED] ■ 

In statistical physics, a coarse-grained Hamiltonian determines the relative 
probabilities of different field configurations for near-equilibrium systems. Ex¬ 
cursions from equilibrium leading to dynamic fluctuations follow by perturbing 
the system with noise. The dynamic response of the system is determined by a 
stochastic partial differential equation , also known as Langevin equation , which 
can be derived in the framework of linear response theory UB|. The correlation 
functions are shown to obey the fluctuation-dissipation theorem , which connects 
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them to respective susceptibility functions [12] ■ In physics, this approach has 
been applied to study dynamic critical phenomena EH HZ]. 

The physical insight provided by the the fluctuation-dissipation theorem is 
that the fluctuations of the unperturbed system contain information about the 
response to external perturbations that drive the system away from equilibrium. 
This formalism, to our knowledge, has not been applied to generate space- 
time covariance functions. Our goal in this paper is to show that effective 
Hamiltonians with local interactions, which are successfully used in spatial data 
analysis, can be extended to space-time random fields. This contribution is part 
of an ongoing effort to transfer ideas from statistical physics to space-time data 
analysis [251 |25] [TU] |25] 1251157] ■ Using the theory of linear response, we show 
that the local interaction framework leads to explicit forms for new space-time 
covariance functions. 

The remainder of the paper is structured as follows: In Section [2] we present 
mathematical background on Gaussian random fields. Section [3] reviews local 
interaction random fields and demonstrates that the so-called Spartan spatial 
random field model is derived from the principle of maximum entropy. Section [4] 
briefly presents linear response theory focusing on the calculation of correlation 
functions. In Section [5] we use linear response theory to derive equations of 
motion for space-time SSRF covariances, and we obtain an explicit equation 
for the SSRF spectral density. Section [6] gives a standard derivation of space- 
time covariances from Langevin equations. This leads to a general equation 
for the space-time covariance function of fields driven by colored noise and 
recaptures the results of the previous Section in the SSRF case with Gaussian 
white noise. In Section [7] we present explicit equations for SSRF space-time 
covariance functions. Finally, we present our conclusions in Section [8] 


2 Mathematical Preliminaries 

In the following, we focus on Gaussian space-time random fields. These can 
be used to model more complicated distributions in the Bayesian framework of 
Gaussian random processes [15]. A space-time random field (STRF) X(s, £;w) € 
R where (s, t) £ R d x R and u> £ is defined as a mapping from the probability 
space (f2, A 1 P ) into the space of real numbers so that for each fixed coordinate 
(s, t) , A'(s, t ; u>) is a measurable function of u> [5]. An STRF involves by definition 
many possible states indexed by to [5] [55]. In the following, we drop the state 
index to for simplicity; instead, we use the symbol w for the cyclical frequency 
in Fourier transforms. The states (realizations) of the STRF are denoted by 
the real-valued scalar functions x(s,t). In the following, <b[Af(s,f)] denotes a 
functional of the STRF X(s,t) that takes unique values for each realization 
x(s,t). 

The expectation over the ensemble of STRF states is denoted by the angle 
brackets (•). Hence, the covariance function is given by 

U x (s, t; s', t') = (X(s',t')X(s,t))-(X(s',t')) (X(s,t)). (1) 
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An STRF is called statistically stationary in the weak sense (or simply sta¬ 
tionary for brevity), if its expectation is constant and its covariance function 
depends purely on the spatial and temporal lags, r and r, respectively. For 
simplicity, since we aim to calculate space-time covariance functions, we assume 
a zero-mean STRF, i.e. (X(s,t)) = 0. The STRF variance will be denoted 
by cr 2 . Furthermore, a stationary random field is called statistically isotropic 
if its covariance function depends only on the Euclidean distance but not the 
direction vector. In the following we focus on statistically isotropic STRFs. 

In the spectral domain, we use the wavevector k to denote the spatial fre¬ 
quency and w to denote the cyclic frequency with respect to time. For a given 
function x(s, t) with space and time dependence, we will use x (k, t) to denote the 
spatial (i.e., with respect to the space variable) Fourier transform and i(k, w) 
the full Fourier transform with respect to both space and time. The pairs of the 
direct and inverse Fourier transforms, respectively, are defined as follows: 

x(k, t) = ,F a [a;](k, t) = f cfee _3kr x(s,t), (2a) 

Jm. d 

x(s,t) = J' 1 7 1 [i](s,t) = — dke jW r x(k,t), (2b) 

( Z7T ) J R d 

x(k, uj) = T B ,t[x\{k, w) = f ds j dte~ J ^' r+ult ^ x(s,t), (2c) 

JR d J — oo 

x(s,t)=J'ifl > [x](s > t)= * f dk f (W (kr+aJt) i(k,u;). (2d) 

[ZTT) t j Rd j _ qo 


3 Local Interaction Random Fields 


Time-independent random fields represent the long-time, i.e., equilibrium, con¬ 
figurations of dynamic variables. They also model “quenched randomness” 
which characterizes the structure of geological media and is a key factor for 
subsurface physical processes |5]. Static random fields can be defined in terms 
of a “pseudo-energy” functional %[a;(s)] which assigns different “energy” levels, 
and subsequently different probabilities, to different configurations x(s). In the 
following, we will assume that the energy functional takes real, non-negative 
values. The most probable configuration of the spatial random field minimizes 
n[x(s)}. 

The joint probability density function (pdf) of the equilibrium spatial random 
field X(s;w) is determined from the functional %[x(s)]. The joint pdf for the 
realization (state) x(s) is proportional to exp[—%[x(s)], i.e., 

/ eq [x{s)] = ±e- H W>\ 


where the partition function Z is given by the following functional inte¬ 
gral [m ng 

3 -W[x(s)] 


Z = J Vx(s) e 
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The functional integral is the continuum limit of the discrete representation 
T>x(s) = limn^oo nr =1 dxi, where the vector (x±,. .., x n ) T represents a dis¬ 
cretization of the continuum state x(s) at n points. Note that whereas Z may 
diverge for n — > oo, the SRF statistical moments are nonetheless well defined. 
For example, the covariance is given by the following functional integral 

(X( Sl )X(s 2 )) = | J Vx(s)x( Sl )x(s 2 )e- n ^. (3) 

In statistical mechanics, the 'H[x(s)] for a given system is obtained from 
kinetic and potential energy terms that reflect the motion and interactions of 
microscopic constituents. In the case of macroscopic random fields X(s ,i), the 
dynamics may not be fully known or solvable. Then, %[x(s)] represents an 
effective functional that incorporates fictitious interactions defined by means of 
derivatives of the field realizations. Hence, %[ x(s)] may represent a sufficiently 
flexible heuristic functional and not a first-principles Hamiltonian. 

3.1 Maximum Entropy Formulation of Spartan Spatial Ran¬ 
dom Fields 

Based on the concept of the effective energy functional, the Spartan spatial 
random field model was proposed [Ml EH EH E3. This model is an extension 
of the classical Gaussian field theory [35] which incorporates a square curvature 
term. The respective energy functional is given by 

H ig c[x(s)\ = J^ds {[z(s)] 2 +7 ?i£ 2 [Vz(s)] 2 +£ 4 [V 2 z(s)] 2 j .(4) 

In the above, 770 represents a scale coefficient, 771 a rigidity coefficient, and 
£ a characteristic length; V stands for the gradient and V 2 for the Laplace 
differential operator. The energy functional TLf gc [x(s)] generates Gaussian, zero- 
mean, stationary random fields. As in classical field theory, a frequency cutoff 
is required for differentiability of the field states [55]. In the absence of such a 
cutoff, the high-frequency contributions lead to random fields that are mean- 
square continuous but non-differentiable in d = 2, 3 [24] ■ Explicit expressions for 
the covariance functions in d = 1,2,3 at the limit of infinite cutoff are derived 
in EX E3 ■ The use of © in spatial data analysis requires the estimation of the 
coefficients 770 , 771 and £ from the data. 

The pdf with energy functional Q can be derived using the principle of 
maximum entropy |32l . The main idea underlying maximum entropy is that the 
pdf is determined by maximizing the entropy of the system given constraints 
imposed by the data. Furthermore, it is assumed that the data constraints are 
equal to the respective expectations over the pdf. 
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Let us define the integrals Sj[a:(s), i = 0,1, 2, by 


S 0 [z(s)] :=j 

( dsx 2 (s), 

(5a) 

s i Ns)] -=j 

( ds [Vx(s)] 2 , 

(5b) 

S 2 [z(s)] :=j 

f ds [V 2 x(s)] 2 , 

(5c) 


In addition, we assume that the expectations ( Sj[a;(s)] ) can be estimated from 
the data. If we possess only one realization, e.g., x*(s), the above assump¬ 
tion requires that the system is ergodic with respect to the constraints, i.e., 
<Si[z(s)]> = Si[x*(s)}. Then, the maximum entropy pdf conditional on these 
constraints is given by 

f [x(s)] = e -/i-AoS , 0 [x(s)]-AiSi[a:(s)]-A 2 ,S 2 [x(s)]. (g) 


The constant /j, normalizes the pdf, i.e., Z = exp (/i), whereas the constants A j, 
i = 0,1, 2 are in principle obtained by solving the system of the three equations 
<Si[s(s)]> = Sj[o;(s)]. Note that Q is equivalent to Q if Ao = l/2r) 0 t; d , Ai = 
? 7iC 2_d /2??o, and A 2 = ^~ d /2r} 0 . 

For data sampled on discrete supports, the S', [a;(s)] are replaced by dis¬ 
cretized estimators S)(x*), that involve the sample values (x*) T = (a:*,... x^)- 
For example, on a regular grid the derivatives are replaced by respective finite 
differences which thus implement the high-frequency cutoff (2Ti| . On irregular 
grids, estimators can be constructed using kernel functions MM- Then, the 
estimation of the model parameters is reduced to the solution of the following 
moment equations 


S 0 (x*) 

Si(x*) 

£ 2 (x*) 



dxN S 0 (x) e -^- A o«o(x)-A 1 S 1 (x)-A 2 5 2 (x) ; 
dx N Si(x) e -^- A o®o( x )- A i^(x)-A 2 5 2 (x) ; 



dx N S 2 (x) e -M->‘oSo(x)-A 1 S 1 (x)-A 2 S 2 (x)_ 


(7a) 

(7b) 

(7c) 


4 Dynamic Correlations 

If a perturbation drives the system away from the equilibrium and the deviation 
is small, the equilibrium fluctuations determine the non-equilibrium response. 
This idea is exploited in statistical mechanics by means of the linear response 
theory [251 . 


4.1 Susceptibility Function 

The theory of linear response focuses on small perturbations of the field X(s, t) 
around the equilibrium state which are caused by an external field h(s,f). If we 
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define by (-) 0 the expectation over the equilibrium distribution and by (•) the 
non-equilibrium expectation, the response of the system to the external field 
can be expressed as follows in terms of the susceptibility function x(s — s'; t — t') 
and the following convolution equation 


(X(s,t)) = (. X(s,t)) 0 



ds' y(s — s'; t — t') h( s', t'). 


( 8 ) 


For large perturbations, the response should also include nonlinear terms. 

Since the difference (X(s,t)) — (X(s,t)} 0 = (X'(s,t)) corresponds to the 
expectation of the non-equilibrium fluctuation, Eq. ([8]) implies that the suscep¬ 
tibility is given by the following functional derivative 


%(s — s' ,t — t') 


lim 

^(s'.tq-s-o 


S(X'(8,t)) 
Sh( s', t') 


(9) 


In the above, 5(-)/5h(s',t') represents the functional derivative, which is defined 
by means of 5h(s, = 5(s — s') 8{t — t’). The limit h(s', t') —► 0 is taken 

in order to ensure that nonlinear terms vanish. 

Let A' 7 (k, w) and h(k, u>) represent the Fourier transforms of the STRF fluc¬ 
tuation and the external field which are given by <©• The Fourier transform of 
the susceptibility function is then given by the following limit 


6{X'{ k,w) ) 

X(k,a;) = ]im . (10) 

ft.—>-o o/i(k, u!) 

The susceptibility function is used in physical applications to describe the 
system’s response to measurable external fields (e.g., electric or magnetic fields). 
In the analysis of spatial data, however, the external field does not necessarily 
represent a physical reality; hence, the notion of susceptibility has not found 
wide applicability. Nevertheless, we suggest that there are potential applications 
of the susceptibility concept, since the analysis of spatial fields often involves 
auxiliary variables [52], which essentially represent external fields. For example, 
the elevation is an auxiliary variable that has significant impact on rainfall; thus, 
it makes sense to consider the susceptibility of the rainfall field to the elevation 
of a particular site. 


4.2 Langevin Equations 

Within the framework of linear response theory, the STRF dynamics is derived 
from the equilibrium energy functional and the noise held that drives deviations 
from the equilibrium by means of the following Langevin equation mm 


dX(s, t ) 

Ft 


-r 


6W.lt e ( s )] 
6x( s) 


tc(s)=X ( s,t) 


C(s,t), 


( 11 ) 
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where T is a diffusion coefficient, <5[.]/<Jx(s) denotes the functional derivative with 
respect to the field state, and £(s, t) is the noise field. The latter is typically 
Gaussian white noise with (£(s, t )) = 0 and variance equal to D 1 i.e., 


(C(s, t) C(s', t')) = D8{ s - s') 8(t - t'). 


( 12 ) 


Equation (11) links the rate of change of the field state x(s, t) to an equilibrium¬ 
restoring velocity that depends on 'H[x(s)] and a random velocity given by the 
noise term. For example, if "H[x(s)] is given by the Ginzburg-Landau effective ac¬ 
tion, then (11) is known as the time-dependent Ginzburg-Landau equation |35l 


p. 192]. In the following, we will denote the restoring velocity for a given random 
field state by 

<13) 

and C[X(s,f)] will denote the respective functional of the random field X(s,t). 


4.3 Fokker-Planck Equation 

The pdf of the STRF X(s,t) that is governed by the Langevin equation © is 
the solution of the following Fokker-Planck equation m 


df [cc(s, t)] 
dt 


J Rd . ds Sx{s,t) . 


-C[x(s',f)] f x [x(s, t)] 


D Sf x [x(s,f)j 
2 Sx(s',t) 


The equilibrium (time-independent) pdf is the asymptotic limit (as t 
the solution of the above Fokker-Planck equation: 


fe q [a<s)] oc exp <^ - 


2TH [x(s)] 
D 


(14) 
00 ) of 

(15) 


The Fokker-Planck equation may not admit an explicit solution; however, 
for many applications in spatial data analysis it is sufficient to know the co- 
variance of the random field, since the fluctuations often follow the Gaussian 
law, whereas in other cases the data can be transformed by means of nonlinear 
transformations to approximately fit the Gaussian law [52] . 


4.4 Equation of Motion for the Covariance 

To obtain the equation of motion (EOM) for the covariance function, we follow 
the approach described in [121 PP- 120-121]. First, we assume that t > t! 
without loss of generality. We use the covariance definition 0 - we replace the 
time derivative of X(s,t) with and then replace x(s, t) with the random 
fields X(s,t). These steps lead to the following equation 


<9C x (s, t; s', t') 
dt 


([E[A'(s, f)] + C(s, t; w)] X(s', t ')). 
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(16) 












The right hand side of the above equation involves the term ( ( (s, i; u) X (s', t')). 
The principle of causality, however, implies that the noise at time t can not in¬ 
fluence the field at the earlier time t'. Hence, noise-field cross correlation is 
dropped, and we obtain the following EOM for the covariance: 

- g~ t - = {V[X(s,t)\X{s .t)). (17) 

Let us now assume that t < t' then since we can not reason that causality 
causes the noise-field cross correlation to vanish the EOM for the covariance 
will be given by: 


9C' x (s,t;s , ,t / ) 

dt 


= ( V[X(s, t)} X(ef, t')) + (C(s, i; w) X{sf, t'-, to)). (18) 


Moreover from (17) the partial derivative of the covariance with respect to 
t' will be given by: 


dC(s, 

dt' 


= (V[X(s , ,t , )]X(s,t)}. 


(19) 


Next, we subtract each side of the equation for <9C x (s, t; s', t')/dt' from the 
respective side of (18) for 9C x (s, t; s', t')/dt. At this point we use the stationarity 
property, i.e., that the covariance is a function only of the lag t—t 1 , which implies 
dC x (s,t;s',t')/dt' = —dC x (s,t-,s',t')/dt. These operations lead to the following 
equation for the covariance rate 

0 <9C x (s, t; s', t') ^ V/J 


dt 


-(y[X(s , ,t , )]^(s,t;cu)) 
+ (C(M;w) X(s\t';u>)). 


( 20 ) 


Using time and space translation invariance, the first two averages on the 
left hand side cancel each other. Combining Equations ©, © and ( [20] ) we 
obtain 


dC{ s, t; s', t') 
dt 


= sign(t - t’) ( V[X(s’, t')] X(s,t )). 


( 21 ) 


4.5 The Fluctuation-Dissipation theorem 

The fluctuation-dissipation theorem links the covariance and the susceptibility 
functions [23]• For t > t! from (20) it is straightforward that: 

dCx{s,t;s’,t') 


dt 


= {((s',t';u)X(s,t;cv)). 


( 22 ) 
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The term in the right hand side of can be evaluated by means of the Furutsu- 
Novikov theorem [TBJ 03]. The latter states that if £(s, t; u>) is a Gaussian process 
and </>[X(s,i)] is a function of X(s,t) then the following identity holds 


(<f>[X(s,t)]{(s',t';u>)) = / ds" 


dt" 


S<p[X(s ,t)} 


(C(s',t»C(sV») 

(23) 


In the present case, 0[X(s, t)] = X(s,t). We use the definition of the sus¬ 
ceptibility function ( |9[ as a response to the noise field, and we take into account 
the noise covariance ( |12[ ). Then, the time derivative of the covariance function 
is related to the susceptibility function via 


0 <9C x (s,t; s',t') 
dt 


D x(s, t, s', t'). 


(24) 


Furthermore, for statistically stationary and homogeneous conditions, the 
above equation becomes 


<9C X ( r,r) 
“ dr 


Dx(r,r). 


(25) 


5 SSRF-based Space-Time Covariance Functions 

Below we derive a partial differential equation for the space-time covariance 
function based on the equation of motion We consider a slightly modified 
form of the Spartan energy functional Q, in which the curvature term is mul¬ 
tiplied by the coefficient fi > 0. In the following, we assume an infinite spectral 
cutoff for the spatial frequency (wavenumber). 

We assume that T = D /2 in line with the relaxation model A of Hohenberg 
and Halperin [Eq.(4.2), p. 445] in [23]. It can then be shown using the properties 
of functional derivatives and ignoring boundary effects [26| that the Spartan 
energy functional is associated with the following restoring velocity 

V[X(s,t)\ = --£- (1- ??1 £ 2 V 2 + ^ 4 V 4 ) X(s,t). (26) 

2 £“ Vo 

In light of the above restoring velocity and ( fl7| ), the stationary covariance 
function is given by the solution of the following equation of motion 

dC f T ,T) = -bsign(r) (1 -VifVl + m£ 4 V 4 ) C x (r,r), (27) 
where D = is the combined diffusion coefficient. 

2£ d i)o 

We define by C(k, r) the Fourier transform of the covariance function over 
the spatial lag vector. This is given by the following d-dimensional integral 

C( k, r) = J- r [C](k, r) = f dr e _jk r C( r, r), (28) 

J R d 
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The function (7(k, r) is the spectral density of the covariance at time lag r. 
Expressing the equation of motion (27) in Fourier space, we obtain the following 
spectral counterpart 


dC ^ T ^ = --Psign(r) (1+rnk 2 ^ 2 + /rfc 4 £ 4 ) C(k,r), (29) 

where k = ||k||. The first-order ordinary differential equation ( [29] ) admits the 
following exponentially damped time evolution of the Fourier modes 

<5(k, r) = C{ k, 0) |t| fc 4 « 4 ) ; (30) 

where C'(k.O) is the covariance spectral density at zero time lag. 

Let us assume statistical isotropy of the spectral density at zero time lag, 
i.e., (7(k, 0) = C{k, 0). Then, based on the spectral representation , the space- 
time covariance function is given in real space by the following one-dimensional 
integral [55] 


A-d/2 „-D |r| /-oo 


r oo 

C{ r,T) = ' {2 ^ d/2 J o &C(M)^V 1 We- DH(,1 ^ 2+ ' ltr ), 

(31) 


where Jd./ 2 -i( x ) is the Bessel function of the first kind and of order d/2 — 1 . 
Given that the covariance function is stationary, the expression (31) also holds 
for t —> oo. At zero lag, therefore, C(r, 0) should be identical to the spatial 
covariance function that corresponds to the static energy functional J§- By 
setting r = 0 in (31) it follows that the above constraint is satisfied if C(k, 0) is 


C(k, 0) = 


riot a 


1 + m(k€) 2 +ia{k& 4 


(32) 

slightly modified by 


which represents the SSRF spectral density given by 
the presence of /i. According to Bochner’s theorem [3], (32) is a permissible 


spectral density if (i) C(k, 0) > 0 and (ii) f dkC(k,0) < oo. These conditions 
are trivially satisfied if rji > 0, p > 0 and d < 3. 

The space-time spectral density S(k,w) is defined as the Fourier transform 
of the space-time covariance function in the spatial-temporal frequency domain, 
i.e: 


S(k,w) =^ P , T [C](k,w)= /" dr f dre~^ r+ulT) C(r,T). (33) 

J —oo J 

The spectral density 5(k, u>) can also be expressed by means of the spectral 
density of the spatial fluctuations at time lag r which was obtained in (30) as 
follows: 


/ oo r oo 

dre~ 3U T C(k, r) = C(k, 0) / dre- J “ T e~ D |r| (i+mk a e+pk « 4 )_ 
-oo J — OO 

( 34 ) 
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Figure 1: Space-time spectral densities for d=l (a) and d=3 (b) and for 771 
1, r)o = 1, £ = 3, D = 0.5, m = 1 
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Using the isotropic C'(k, 0) given by (32) and taking into account that the 
integrand in (34) is even, 5(k, w) is expressed in terms of the following cosine 
Fourier transform 


S(k,w) 


_ 2? 7o£_ 

1 + ?7ifc 2 £ 2 + n fc 4 £ 4 



dr cos (uit) 


e -D\r\(i+ m k 2 e+^kU 4 ) 


(35) 


For rji > 0 and D > 0 the following analytic expression for the the SSRF 
space-time spectrum is obtained using the identity [Eq.(1.4.1) p.14] in [41] 
Eq.(1.4.1) p. 14] 


S(k, w) = 


2??o £D 


D 2 ( 1 + rjik 2 ^ 2 + n fc 4 £ 4 ) 2 + uj 2 


(36) 


The above is a permissible spectral density function of the form S(k,ui) = 
a/(II(k) + ui 2 ), where II(k) is an eighth order polynomial. This model is more 
flexible than the one described by the second order Matern spectral density, 
commonly used in geostatistics, in which the respective polynomial II(k) is 
given by II(k) = (k 2 + cj ) 2 ) 4 [33]. In Figure 0 the spectral density ( [36] ) is 
plotted for d = 1 and d = 3 for fj, = 1. The spectral density (36) could by used 


in meteorological modelling as space-time extension of purely spatial spectral 
densities derived by polynomials of the diffusion operator m- 


6 Covariance Function from Langevin equation 


In this section we obtain space-time covariances from the equation that governs 
the motion of the STRF realizations. Based on (11) and ([26]) the following 


Langevin equation governs the time evolution of the state x(s,t) 


dx(s,t) 

dt 


= V[x(s,t)] + ((s,t), 


(37) 


where V[x(s, t)] is the restoring “velocity” as defined by (13), and ((s,t) is the 
stochastic “velocity”. We assume a bilinear equilibrium energy functional of the 
form 


H [x(s)[ = f ds f ds' x(s') C(s' — s) x(s), 

Jm d Js. d 

where C(s' — s) is the inverse spatial covariance operator defined by means of 
the convolution integral [54] 


f dsi £(s — Si) C(si 

J R d 


s') = S(s — s'). 


13 











Note that £( s' — s), M d i— > ffi. is a real-valued, and symmetric function, i.e., 
£(s' — s) = £(s — s'). The restoring velocity is given by 


V[x(s,t)] 


dsi £(s — Si) x(s i,t). 


The Langevin equation corresponding to (37) for the spatial Fourier modes of 
the STRF is 

3x(k, t) 


dt 


= -£(k)i(k,f) +C(k,t), 


where £(k) £ K is the inverse spatial covariance operator in Fourier space. 
If £(•) is a linear combination of even-order derivative operators, then £(k) 
respectively is a polynomial in k 2 . 

The temporal evolution of the Fourier modes is given by the solution of the 
above ordinary differential equation, i.e., by 


s(k,t) = x(k,0)e- £(k)t + [ dt 1 e" £(k) C(k, £'). 

Jo 


(38) 


The STRF spatial spectral density is given by C( k, t: k', if) = (x(k, t) x(k', t!) ) c . 
It follows from the above that the temporal evolution of (7(k, t; k', t') is deter¬ 
mined by the integral equation 

(®(k, t) i(k', H) ) c = e - £ ^ 4 e" f < k '> 4 ' (®(k, 0) s(k', 0) ) c (39) 

+ J dti j dt 2 e~ Z W e -^( k ') (*'-*») ^ C(k, tr) C(k 7 , * 2 ) ), 

where (x(k, 0) x(k', 0)) is the initial condition at zero time lag. 

If we assume a correlated stochastic “velocity” with spatial spectral density 


( C(k, h) C(k', t 2 ) ) = 2 D (2n) d <5(k + k') c c (k) S(h - t 2 ), 


we obtain the following integral equation for the temporal evolution of C(k, t ; k', t r ) 
C , (k,i;k',i') = e-[ f(k)t+f(k ' )t ' 1 C(k,0;k',0) (40) 

+ (2tt ) d D / dtie“[ /:(k) ( t - tl)+£ ( k ' ) ( t, “ tl) U(k + k , )c c (k). 


Setting t = t' in (40) we obtain 


C(k,t-, k',£) 


e -[ £ (k)+ £ (k')]‘c>( k ,0 ; k',0) 


+ (27r) d D 


5(k + k') c^(k) 
£(k) + £(k') 


\ _ e -(£(k)+£(k'))t 


(41) 
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Due to stationary it holds that C(k,t;k' ,t) = C(k, £; k', t') for all t and t'. 
Hence, we obtain the following relation for the spatial spectral density at zero 
time lag 


C(k, 0; k', 0) = lim C(k,f;k',f) 

t—t OO 


(2 Tr) d D 


<5(k + k') d( (k) 
£(k) + £(k') 


(42) 


Combining (401 and (421, the temporal evolution of the spatial spectral density 
is given by the following exponential relation 


C(k,t;k',t') = C( k j O;k',O)e- [0(t -‘ ,)£(k) ( t -*') +e ( t '- t ) £(k ')( t '-‘)l. (43) 


Finally, in light of (43) and by invoking the reflection symmetry £(k) = £(—k), 


we obtain the following equation for the evolution of the space-time covariance 


C( s, t; s', t') = 


dk 


dk' 


'R d (2n) Js. d (27t)° 


e J ( ks + k ' s ')c'( k ,t ;k , ,£) (44) 


D 


f dk 

m d (2n) a 


e jk (s-s') e -C(k) 


C£( k ) 


Equation (44) holds in general for colored noise and for general £(fc). For 


the SSRF case, the inverse spatial covariance operator £(fc) is given by [27 
A k ) = [! + Vik 2 ^ 2 + M^ 4 ] ■ 

Assuming that the stochastic “velocity” is Gaussian white noise, i.e., c£(k) = 1, 
the SSRF covariance C(r, r) is given by 


C(r;r) = 


dk 


0 J kr . 




(2n) d 1 + di* 2 £ 2 + d 


e -^( 1 +J ' lfc 2 « 2 +^ fc 4 « 4 )M. ( 45 ) 


Based on the above and the spectral representation of isotropic random fields 
which transforms the multidimensional integral into a one-dimensional inte¬ 


gral [55], we recover (31) derived above. 


7 Explicit Space-time Covariance Functions 

Equation © is the general result for the space-time covariance function ob¬ 
tained in the framework of linear response theory from the equilibrium local- 
interaction energy functional Q. We use the following integral identity for the 
spectral density 




1 + 771 k 2 £ 2 + ii k 4 £ 4 J o 


/>00 

/ dne~ K [ 1 + , n( fe £) 2 +M fe £) 4 ]/( 7 ?o£‘ i )_ 

Jo 


(46) 
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The above is based on 1/A = due kA , for A > 0. 

The above leads to a space-time covariance function which is given by 




r l-d/2 e -D |t| 
(27t)' 1 / 2 




dkk d ' 2 J d/ 2 -i{kr)e~ A ^ k \ 


(47a) 


A(n,k)=D\r\ (r)ik 2 £ 2 +/J.k 4 £ 4 ) +n [l + ??i(fc£) 2 + M(*0 4 ] /(Vo £ d )- ( 4 7b) 

We can thus express the space-time covariance function in terms of the following 
double integral 

l-d/2 -D\t\ poo 

C ( r > r ) = {2n)d /2 J o dnQ(K), (48) 

/*oo 

Q(k)=/ dkkWja/i^We-^*’®. (49) 

Jo 


7.1 Covariance for the zero-curvature model 


We investigate the space-time covariance (47) with p = 0. In this case the 
function A(n,k) is given by 


A («, k ) = + (D\t\ m f , 2 + 


KVlZ 2 

vof d 


k 2 . 


(50) 


Then, the spectral integral (49) can be analytically performed using [Eq.(11.4.28), 
p. 486] |2] 


» T - ' ( 

dk e~° 2 fc2 k^~ l Jv(bk) = 1 2 ’ 


2a ^ T(i/ +1) \2a 


x- 1 i F l 


v + n b 2 

2 > ^ + 1; — 4^2 


(51) 

In the above, the iFi(ai, 02 ; z ) —also known as M(a\, a 2 ',z )— is the conflu¬ 
ent hypergeometric function, defined by [Eq. (13.1.2), p. 504] [2] as follows 


M(a 1 ,a 2 ;z) = ^ 


n —0 


(ai ) (n h n 
(02)^ n\ ’ 


where (ai)(”) is the rising factorial, 

(ai) (n) = ai(ai + 1)... (aq + n - 1), (a)[ 0) = 1. 

We apply the above definition with a 2 = D |r| rft £ 2 + k r)i f 2 /(r] 0 £ d ), v = 
d/2 — 1, p, = d/2 + 1, and b = r to the integral in (51). Based on these 
replacements, (49) and (50), we evaluate the function Q(k) as follows 


Q(k) 


e / r \ d/ 2 —i /d d — r 2 \ 

2a d / 2+1 \2a) 1 1 ^2’ 2’ 4oV ' 


(52) 
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Further, using the definitions for the hypergeonretric function, it follows for 
a\ = a ,2 = d/2 that M{d/2,d/2-,z) = exp(^). In light of this identity, and the 
function Q{k) is given by 


r d/ 2-1 




(53) 


Hence, the expression (53) for Q(n) and (48) lead to the following univariate 
integral for the space-time covariance function 


C(t,t) = 


M e ±(f>i+f>2K) 


-k/3 0 


(47r) d / 2 


dn 


{Pi + P 21 P) 


d/2 


where 


Po = Pi=D \t\ (? 7 iC 2 ) , P 2 = 


m/r 

Vo C 


d ’ 


(54a) 


(54b) 


Equations (31) and (54) both involve univariate integrals. The latter, how¬ 
ever, involves a non-negative, decreasing function of k, whereas the former in¬ 
volves the integration of the oscillating Bessel function Jd/ 2 -i{kr) which is more 
difficult to evaluate numerically, especially for large lags. 


7.1.1 The case of zero spatial lag 


For r = 0 the covariance function (54) is given by the integral 

OO p-K-Pt) 


C(0,t) = 


e~ D l r l 
(47r/3 2 ) d / 2 


dn - 


M d/2 ' 
02 J 


(55) 


The integral can be performed using integral tables, and more precisely 
[Eq. (3.362.2), p. 362] [19] for d = 1, [Eq. (3.352.4), p. 358] £H] for d = 2 and 
[Eq. (3.369), p. 363] [19] for d = 3 which lead to the following expressions 


C{0,t) = { 


-S^Ei(-D 

2rioe~ D T _2 m__ 

?r( 4 »Ji) 3/2 



(56) 


In the above, Ei(:r) is the exponential integral function defined by Ei(:r) = 
Jdte*/t and Erfc(a;) is the complementary error function defined by Erfc(:r) = 

c?te _t . The expansion of Ei(x) around zero is given by 


Ei(z) = 7 + In M + Tf 

fe =1 


where 7 is the Euler-Mascheroni constant. In light of the above and (56) it 
follows that the time covariance has a singularity at r = 0 for d = 2,3. 
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7.1.2 The case of zero time lag 


For r = 0 it follows that fa 
the following integral 


0 and the covariance function (54) is given by 


C(r,0) 


1 

(47T/3 2 ) d / 2 



0 4 02 


:-«A) 


i d / 2 


fa 

(<infa ) d / 2 fa Jo V d/2 


(57) 

where y = fan and /i = r^/ fa/ fa = This integral can be evaluated 

using the identity [3.471.12, p. 385] [TH] which leads to 


<■> 

where K'd/ 2 -i(') is the modified Bessel function of the second kind of order 
d/2 — 1. For d = 1 the covariance function depends on K_ l / 2 (-) = K 1 / 2 (-), for 
d = 2 it depends on K 0 {-), whereas for d = 3 the function K 1 / 2 {-) is involved. 
In both d = 1 and d = 3 the behavior of the Bessel function near zero is 


^±i/2(z) 


^(f) 

— ln(a;/2) 


1/2 

5 

^ 7) 


d= 1,3 
d = 2. 


Hence, in light of (58) the spatial covariance has a singularity at r = 0 in 
d = 2,3, in agreement with the limit of the time covariance as well. This is 
not surprising since the model described above is identical to the Gaussian field 
theory, which is known to have a singular behavior at the origin if a finite cutoff 
is not used [12J p. 227]. 


7.1.3 Explicit covariances for general lags 


For d = 1,3, /i = 0 it follows from (31) that the corresponding covariances are 
given by the equations 


C(r ,t) 


e DN? lo r dk cos ( fcr ) c -D\r\(mk 2 C 2 ) 

* Vi £ Jo ^72 + fc 2 


(59) 


and 


C( r,r) 


e-^ |T| 7? 0 g 
27 t 2 ? 7 i r 


rdk k r {kr \ e- 6 

Jo 


(60) 


The above integrals can be evaluated analytically for rft > 0. Applying 
the identities [ Eq.(1.4.15), p. 15] [?T] and [ Eq.(2.4.26) p.74] UT|, the following 
explicit formulas are derived for the one- and three-dimensional covariances 
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Figure 2: Space-time covariances for (a) d=l and (b) d=3. In both cases 771 = 
1 ,770 = 1, £ = 3, D = 1. In (a) the minimum lag — in both space and time— is 
equal to the floating point relative accuracy, whereas in (b) it is equal to 0.1 to 
avoid the singularity at ( 0 , 0 ). 




(a) (b) 


Figure 3: Space-time covariances for d = 1 ,770 = 1, £ = 3, D = 1. (a): Covari¬ 
ances for zero space lag and different 771 are shown, (b): Covariances for r = 3 
and different rji. 
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(a) 


(b) 


Figure 4: Space-time covariances for d = 1, 770 = 1, £ = 3, D = 1. (a) Spatial 
dependence at r = 0. (b) Spatial dependence at r = 3. 


C{ r;r) = 


?7o 

4 ^ 


e v^tt« Erfc \/D Irl — 




(61) 


+ e' flri Erfc \ D\t\ + 


2\ Din |r|£ 


and 


C( r ; t) = 


n o£ 

r 


e Erfc I 


-e^i s Erfc 



(62) 


2\/E ) ?7i |t|£, 


2ril?7 ? i |r|£ 


The covariance function (61) is equivalent to a model derived by Heine start¬ 


ing from the general form of a second-order stochastic partial differential equa¬ 
tion ED [S3]- Our result, however, is expressed in terms of the parameters 
Vo-i Vit D that can be easily identified with properties of the STRF and in¬ 
clude the SSRF parameters of spatial random fields. In particular, 770 is an 
overall scale coefficient that determines the amplitude of the fluctuations and 
has units [A '] 2 ,771 is a dimensionless rigidity coefficient that suppresses large gra¬ 
dients, £ is a characteristic length, and D is a characteristic inverse time. We 
also derive the covariance function (62) which is valid in three spatial dimen¬ 
sions. Note that the presence of 771 in addition to £ implies that the correlation 
length is a function of both parameters [3D]. In both d = 1 and d = 3 for r —► 0 
and r —> 0 , the purely temporal and purely spatial covariance models, respec¬ 


tively (561 and (58), are recovered from (61) and (62). In agreement with the 
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zero-lag analysis, the singularity at r = 0 and r = 0 is present for d = 3 in (62). 


This singularity, which denotes that the STRF has infinite variance, is known 
in Gaussian field theory as well as in statistical modelling [83] . We discuss a 
remedy for this problem in the next section. 

In Figure [2] we plot the covariance functions |6~i~[ ) and (62) for d = 1 and 
d = 3, respectively. For d = 3 the minimum lags used are equal to 0.1 to 
avoid the zero-lag singularity. Figures [3] and [4] demonstrate the dependence 


of the covariance function (61) in d = 1 versus the spatial (time) lag under 


constant time (space) lag for different values of the rigidity coefficient r/i. Note 
that at the larger lag distances, r = 3 and t = 3 the covariance function is 
smooth at the origin, i.e., for r = 0 and r = 0, respectively, in contrast with the 
non-differentiable peak at (r = 0,r = 0). 


7.2 Small-/! approximation 


In the following we assume that /r <C 1, so that the curvature term in (47) is 


present but multiplied with a small coefficient. This small perturbation, how¬ 
ever, is sufficient to tame the singularity at the origin (zero lag) of the covariance 
function, effectively by reducing the impact of large-frequency fluctuations. 


We use the spectral representation (47) for the covariance function. The 
exponential e _j4< - K,fc ) in this case is expressed as follows 


g—A(K,k) _ «/ 3 q —u 2 k 2 —fivk 4 


where the coefficients v and u 2 are defined by means of the following relations 


and the constant coefficients given by (54b I 

v = £ 4 (t>|t| + Kfiaj , u 2 = ? 7 i£ 2 (f) |t| + k/3 0 ) . 


(63) 


The exponential component exp(—/t v k 4 ) can be approximated using the Taylor 
expansion truncated at an even order 2 M as follows 


2 M 


exp (-/ xv k 4 ) « Y] 


(-l) m (nvk 4 ) 7 


m—0 

The even truncation order ensure that the approximation is non-negative. The 
divergence of the expansion for large k is controlled by the exponential exp(— u k 2 ) 
for rji > 0. Therefore, the truncated spectral density corresponds to a permissi¬ 
ble covariance function. The resulting expansion of the covariance function ([47 
is given by the following series of integrals 


1—d/2 —ri|r| 2M , 


m—0 


(27r) d / 2 

nOO 

Qm(n) = e-^° / dkk d / 2+4m J d/2 _!(fcr)e 
J 0 


dn v m Q m {n), 


Z 1 z 

—u k 


(64) 

(65) 
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The integral over k in Q m { k) can be evaluated using (511 which leads to 


^ e _K/3 ° (r\ d ! 2- 1 T(d/2 + 2m) 

Qm{K] ~ u 4^ UJ T(d/2) 


, i d —r 2 

1 F\ ( 2?7r + —, —; —— 
111 2’ 2 4« 2 


( 66 ) 


Based on the two equations above, the space-time covariance function be¬ 
comes 




2e" £> l r l 


2 M 

E 


{-fi) m T(d/2 + 2m) 


(4 7r )d/2 ^—' m! 

' ' m=0 


nOO 

Rm{r,r) = dK 
Jo 


e ~K0o v m 


T(d/2) 

d d — r 2 


4m+d li?1 i 2m+ 2’2’4u 2 


Rm(r, r), (67) 

( 68 ) 


In the above, the coefficients t; and it 2 are defined in (63), whereas the coefficients 


Po,Pi ,^2 are defined in (54b). 


7.3 Numerical solutions for general values of fi and rji 


The presence of the parameter /i in the general space-time model (31) eliminates 


the observed singularities at the origin for fj, = 0 . For /j, > 0 and negative 
values of the rigidity coefficient 7 / 1 , i.e. for —2 < rji < 0 the covariance function 
develops oscillations. 


Based on (|3l|, the space-time covariances for fjL ^ 0 in d = 1, 3 are given by 

D\r\( m k 2 f+n.k 4 i 4 )^ 


the following one dimensional integrals 

? 7 o£ f°° j, cos(fc?' 


C(r,r)=- 


dk 


l + 7/i k 2 £ 2 + fi fc 4 £ 4 


C(r, t) = 


e-^ |T| 7? 0 e 3 

27t 2 ?7i 




fc2 sinc(fcr) c _,p| T | (70) 


1 + 7/1 fc 2 £ 2 + M fc 4 C 4 


The above integrals were evaluated in Mathematica using the NIntegrate 
function with high working precision. The spectral integrals were truncated 
using a finite spectral cutoff, i.e. k c = 100. In Figures [5] an d [ 6 ] we plot the 
spatial dependence of the covariance functions (69) and (|70[) for d = 1 and 
d = 3, respectively. Plots correspond to different values of fj and two different 
time lags. In all cases, the correlation length drops with increasing fi. Moreover, 
it is evident in Fig [ 6 ] that in the d = 3 model has finite values at the origin for 

fi > 0 . 

Figure [7] displays how the covariance function (61) varies in d = 1 with r for 
constant r and three negative, values of 7 / 1 . In all cases, the covariances obtained 
exhibit oscillations; such functions find applications in hydrology and in wave 
phenomena E2Z3- In d = 3 for negative values of 7/1 oscillating covariances are 
obtained by (70); in this case, however, the negative hole of the oscillations is 
reduced. 
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(a) (b) 


Figure 5: Space-time covariances for d = l,r] 0 = l,rji = 0.5, £ = 3,D = 1. (a) 
Spatial dependence at r = 0. (b) Spatial dependence at r = 3. 




(a) (b) 


Figure 6: Space-time covariances for d = 3,7^ = 0.5 ,770 = 1,£ = 3,D = 1. (a) 
Spatial dependence at r = 0. (b) Spatial dependence at r = 3. 



(a) (b) 


Figure 7: Oscillating space-time covariances for d = 1, /x = 1, £ = 3, D = 1. (a) 
Spatial dependence at r = 0. (b) Spatial dependence at r = 3. 
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8 Conclusions 


The analysis of spatiotemporal data and statistical mechanics have followed to a 
large extent separate paths, in spite of the conceptual overlaps between the two 
fields. The connection between statistical mechanics and geostatistics has been 
investigated in a series of papers 1231251 I^EBIESIEZI. Herein we demonstrate 
that space-time covariance functions for the analysis of spatiotemporal data 
sets can be obtained using ideas from statistical mechanics, such as held theory, 
maximum entropy, and linear response theory. 

We show that a Gaussian held model with an energy functional to which 
a curvature term is added can be derived from the principle of maximum en¬ 
tropy using suitable data-based constraints. The spatial covariance functions 
obtained from this energy functional incorporate “interactions” between the 
sample locations, and they offer increased hexibility due to a richer paramet¬ 
ric structure [28] [30]. The “equilibrium” spatial model is herein extended to 
the space-time domain by means of the covariance equation of motion or the 
associated Langevin equation for the realizations, both of which are obtained 
within the framework of the relaxation approximation in linear response theory. 
The covariance equation of motion provides explicit formulas for the temporal 
evolution of the spectral density. It is also shown that the associated covariance 
functions in real space involve integrals that can be analytically computed in 
different approximation regimes. The obtained solutions provide new types of 
non-separable, space-time correlation functions without using simplifying but 
unrealistic assumptions. In contrast with purely statistical models, the STRF 
model parameters introduced in this paper correspond to distinct properties of 
the random field. 

For future investigations, flexible representations such as the Karhunen- 
Loeve expansion can be used in combination with the Langevin equations to 
obtain compact expressions for space-time covariance models. The general for¬ 
mulation for construction of space-time covariances based on concepts of statisti¬ 
cal mechanics can be extended to incorporate more general models described by 
non-linear generalized Langevin equations and colored or non-Gaussian driving 

noises [3] 051. 
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